package it.unibo.refolding.pico;

import it.unibo.refolding.common.CommonThings;
import it.unibo.refolding.common.LinearFunction;

import java.awt.geom.Point2D;
import java.util.logging.Level;
import java.util.logging.Logger;

import org.apache.commons.lang3.builder.HashCodeBuilder;
import org.apache.commons.lang3.builder.ReflectionToStringBuilder;
import org.apache.commons.lang3.builder.ToStringStyle;

import weka.core.matrix.Matrix;

import com.google.common.base.Preconditions;

public class LinearRegression {
	
	private static final Logger logger = Logger.getLogger(LinearRegression.class.getSimpleName());
	
	public final LinearFunction f;
	public final double rmsd;

	public LinearRegression(LinearFunction f, double rmsd) {
		this.f = f;
		this.rmsd = rmsd;
	}

	@Override
	public String toString() {
		return ReflectionToStringBuilder.toString(this, ToStringStyle.SHORT_PREFIX_STYLE);
	}

	@Override
	public int hashCode() {
		return new HashCodeBuilder().append(f).append(rmsd).toHashCode();
	}

	@Override
	public boolean equals(Object obj) {
		if (!(obj instanceof LinearRegression)) {
			return false;
		}
		LinearRegression other = (LinearRegression) obj;
		return f.equals(other.f) && rmsd == other.rmsd;
	}

	public static LinearRegression linearRegressionWithPerpendicularOffsets(
			short[] a, int start, int end) {
		double[] x = new double[end - start], y = new double[end - start];
		for (int i = end; --i >= start; ) {
			x[i - start] = i;
			y[i - start] = a[i];
		}
		return fitLineWithPerpendicularOffsets(x, y);
	}
	
	/**
	 * 
	 * @param a
	 * @param start
	 * @param end
	 * @return null when it fails
	 */
	public static LinearRegression linearRegression(short[] a, int start, int end) {
		double[][] data = new double[end - start][2];
		for (int j = end; --j >= start;) {
			data[j - start][0] = j;
			data[j - start][1] = 1;
		}
		Matrix aMatrix = new Matrix(data);
		double[][] b = new double[end - start][1];
		for (int i = end; --i >= start;) {
			b[i - start][0] = a[i];
		}
		Matrix solution;
		try {
		    solution = aMatrix.solve(new Matrix(b));
		} catch (RuntimeException e) {
			if (e.getMessage().contains("singular") || e.getMessage().contains("rank deficient")) {
				logger.log(Level.WARNING, "Could not do linear regression: {0}.", e);
				return null;
			} else {
				throw e;
			}
		}
		logger.log(Level.FINE, "\n{0}", solution);
		assert solution.getRowDimension() == 2 & solution.getColumnDimension() == 1;
		LinearFunction f = new LinearFunction(solution.get(0, 0), solution.get(1, 0));
		Point2D.Double[] points = new Point2D.Double[end - start];
		for (int i = end; --i >= start;) {
			points[i - start] = new Point2D.Double(i, a[i]);
		}
		return new LinearRegression(f, f.rmsd(points));
	}
	
	/**
	 * http://mathworld.wolfram.com/LeastSquaresFittingPerpendicularOffsets.html
	 */
	public static LinearRegression fitLineWithPerpendicularOffsets(double[] x, double[] y) {
		Preconditions.checkArgument(x.length == y.length);
		double xMean = CommonThings.mean(x), yMean = CommonThings.mean(y);
		double B = 0.5 * (sumSquareLessSumMeanSquare(y) - sumSquareLessSumMeanSquare(x)) / (-sumProductLessMeansProduct(x, y));
		if (Double.isInfinite(B)) {
			logger.log(Level.WARNING, "Could not compute the least squares fitting line with perpendicular offsets"
					+ " for arrays of lengths {0} and {1}.", new Object[] {x.length, y.length});
			return null;
		}
		// In contrast to the wolfram page, we use y = a * x + b.
		double a0 = -B - Math.sqrt(B * B + 1), a1 = -B + Math.sqrt(B * B + 1);
		// Perhaps there is a way to know which one is better, or if they are equivalent, but for
		// now just compute the function to be minimised.
		double b0 = yMean - a0 * xMean, b1 = yMean - a1 * xMean;
		double residual0 = residualPerpendicularOffsets(x, y, a0, b0),
		    residual1 = residualPerpendicularOffsets(x, y, a1, b1);
		if (residual0 < residual1) {
			return new LinearRegression(new LinearFunction(a0, b0), residual0);
		} else {
			return new LinearRegression(new LinearFunction(a1, b1), residual1);
		}
	}

	private static double residualPerpendicularOffsets(double[] x, double[] y, double a, double b) {
		double denominator = 1 + b * b;
		double sum = 0;
		for (int i = x.length; --i >= 0; ) {
			double diff = y[i] - (a * x[i] + b);
			sum += diff * diff / denominator;
		}
		return Math.sqrt(sum);
	}

	static double sumProductLessMeansProduct(double[] x, double[] y) {
		Preconditions.checkArgument(x.length == y.length);
		double xMean = CommonThings.mean(x);
		double yMean = CommonThings.mean(y);
		double sum = 0;
		for (int i = x.length; --i >= 0; ) {
			sum += x[i] * y[i] - xMean * yMean;
		}
		return sum;
	}

	static double sumSquareLessSumMeanSquare(double[] x) {
		double a = CommonThings.mean(x);
		double aSquare = a * a;
		double sum = 0;
		for (double xx : x) {
			sum += xx * xx - aSquare;
		}
		return sum;
	}	
}
